Datos abiertos SSA https://www.gob.mx/salud/documentos/datos-abiertos-152127
Diccionarios de datos SSA https://www.gob.mx/salud/documentos/datos-abiertos-152127 y Catálogo Único de Claves de Áreas Geoestadísticas Estatales, Municipales y Localidades de INEGI https://www.inegi.org.mx/app/ageeml/
## Warning: package 'zoo' was built under R version 4.0.3
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: sp
## rgdal: version: 1.5-15, (SVN revision 1045)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: C:/Users/vshal/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: C:/Users/vshal/Documents/R/win-library/4.0/rgdal/proj
## Linking to sp version:1.4-2
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
ruta_estados <- '../datos_shp/Estados.shp'
estados <- rgdal::readOGR(ruta_estados)
## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : Discarded datum International_Terrestrial_Reference_Frame_1992 in CRS
## definition: +proj=lcc +lat_0=12 +lon_0=-102 +lat_1=17.5 +lat_2=29.5 +x_0=2500000
## +y_0=0 +ellps=GRS80 +units=m +no_defs
## OGR data source with driver: ESRI Shapefile
## Source: "D:\GD\Projects_actual\geo_COVID-19\datos_shp\Estados.shp", layer: "Estados"
## with 32 features
## It has 2 fields
ruta_municipios <- '../datos_shp/Municipios.shp'
municipios <- rgdal::readOGR(ruta_municipios)
## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : Discarded datum International_Terrestrial_Reference_Frame_1992 in CRS
## definition: +proj=lcc +lat_0=12 +lon_0=-102 +lat_1=17.5 +lat_2=29.5 +x_0=2500000
## +y_0=0 +ellps=GRS80 +units=m +no_defs
## OGR data source with driver: ESRI Shapefile
## Source: "D:\GD\Projects_actual\geo_COVID-19\datos_shp\Municipios.shp", layer: "Municipios"
## with 2456 features
## It has 4 fields
municipios@data$Clave_Mun_Ent_Texto <- as.numeric(as.character(municipios@data$CVE_MUN)) +
1000 * as.numeric(as.character(municipios@data$CVE_ENT))
municipios@data$Clave_Mun_Ent_Texto <- sprintf("%05d", municipios@data$Clave_Mun_Ent_Texto)
#str(municipios@data)
metadatos_m <- datos_municipio[,1:12]
#str(metadatos_m)
dia_c_m <- datos_municipio[,c(casos_fechas)]
dia_n_m <- datos_municipio[,c(neg_fechas)]
dia_antes_c_m <- cbind(datos_municipio[,c(casos_fechas[-1])],data.frame(inicio = rep(0,dim(datos_municipio)[1])))
dia_antes_n_m <- cbind(datos_municipio[,c(neg_fechas[-1])],data.frame(inicio = rep(0,dim(datos_municipio)[1])))
incremento_m <- (dia_c_m + dia_n_m) - (dia_antes_c_m + dia_antes_n_m)
incremento_c_m <- dia_c_m - dia_antes_c_m
names(incremento_m) <- pruebas_fechas
names(incremento_c_m) <- pruebas_fechas
#str(incremento_m)
incremento_m_a100k <- round(100000 * incremento_m / datos_municipio$Pob_Municipio, 2)
incremento_m[is.na(incremento_m)] <- 0
incremento_c_m[is.na(incremento_c_m)] <- 0
incremento_m_a100k[is.na(incremento_m_a100k)] <- 0
#dim(incremento_m_a100k)
#head(incremento_m_a100k)
semanas_lag <- c(seq(from = 0, to = n, by = 7))
pruebas_breaks <- c(0,10,20,30,40,50,75,100,150,1000000)
pruebas_breaks_labels <- c("0-10","10-20","20-30","30-40","40-50","50-75","75-100","100-150",">150")
pruebas_raw_breaks <- c(0,5,10,20,50,100,200,500,1000,1000000)
pruebas_raw_breaks_labels <- c("1-5","6-10","11-20","21-50","51-100","101-200","201-500","501-1000",">1000")
pruebas_indice_breaks <- c(-1,-0.75,-0.5,-0.25,-0.1,0.1,0.25,0.5,0.75,1)
pruebas_indice_breaks_labels <- c("menos que esperado (-1)","-0.75","-0.5","-0.25","esperado (0)","0.25","0.5","0.75","mas que esperado (1)")
pruebas_positividad_breaks <- c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,1)
pruebas_positividad_breaks_labels <- c("0-10%","10-20%","20-30%","30-40%","40-50%","50-60%","60-70%","70-80%",">80%")
#for (i in c(3,10,30,46)) {
for (i in 1:46) {
selection_cols <- paste0("pruebas",fecha_cadena(as.Date(fecha_inicio) + seq(semanas_lag[i]+1, semanas_lag[i]+8)))
#print(selection_cols)
suma_pruebas_semanal <- apply(incremento_m[,selection_cols], 1, sum)
suma_pruebas_c_semanal <- apply(incremento_c_m[,selection_cols], 1, sum)
total_pruebas_semanal <- sum(suma_pruebas_semanal, na.rm = TRUE)
total_pruebas_c_semanal <- sum(suma_pruebas_c_semanal, na.rm = TRUE)
total_poblacion <- sum(datos_municipio$Pob_Municipio, na.rm = TRUE)
pruebas_positividad_semanal <- suma_pruebas_c_semanal / suma_pruebas_semanal
pruebas_esperados_semanal <- total_pruebas_semanal * datos_municipio$Pob_Municipio / total_poblacion
## indice (similar a NDVI)
pruebas_indice_semanal <- (suma_pruebas_semanal - pruebas_esperados_semanal) / (suma_pruebas_semanal + pruebas_esperados_semanal)
suma_a100k_semanal <- apply(incremento_m_a100k[,selection_cols], 1, sum)
datos_semanales <- data.frame(
Clave_Mun_Ent_Texto = as.character(metadatos_m[,"Clave_Mun_Ent_Texto"]),
Municipio = metadatos_m[,"Nom_Mun"],
Latitud = metadatos_m[,"Lat_Decimal"],
Longitud = metadatos_m[,"Lon_Decimal"],
pruebas_100k = suma_a100k_semanal,
pruebas_raw = suma_pruebas_semanal,
pruebas_esperado = pruebas_esperados_semanal,
pruebas_indice = pruebas_indice_semanal,
pruebas_positividad = pruebas_positividad_semanal
)
#print(str(datos_semanales))
#print(head(datos_semanales, n = 30L))
mis_municipios <- municipios
mis_municipios <- merge(municipios, datos_semanales, by = "Clave_Mun_Ent_Texto", all.x = TRUE)
mis_municipios@data$Class_100k <- as.numeric(cut(mis_municipios@data$pruebas_100k,
breaks = pruebas_breaks))
mis_municipios@data$Color_100k <- brewer.pal(n = 9, name = 'YlGn')[mis_municipios@data$Class_100k]
mis_municipios@data$Class_raw <- as.numeric(cut(mis_municipios@data$pruebas_raw,
breaks = pruebas_raw_breaks))
mis_municipios@data$Color_raw <- brewer.pal(n = 9, name = 'PuBuGn')[mis_municipios@data$Class_raw]
mis_municipios@data$Class_indice <- as.numeric(cut(mis_municipios@data$pruebas_indice,
breaks = pruebas_indice_breaks))
mis_municipios@data$Color_indice <- brewer.pal(n = 9, name = 'PRGn')[mis_municipios@data$Class_indice]
mis_municipios@data$Class_positividad <- as.numeric(cut(mis_municipios@data$pruebas_positividad,
breaks = pruebas_positividad_breaks))
mis_municipios@data$Color_positividad <- brewer.pal(n = 9, name = 'Reds')[mis_municipios@data$Class_positividad]
#print(mis_municipios@data$pruebas_positividad)
#print(head(mis_municipios@data, n = 25L))
#str(mis_municipios@data)
#str(mis_municipios)
## grafica pruebas por 100 mil
plot(mis_municipios,
col = mis_municipios$Color_100k,
axes = TRUE,
border = NA,
#col = "transparent",
main = paste0("Pruebas semanales por 100 mil habitantes ",
as.Date(fecha_inicio) + semanas_lag[i] + 1," - ",
as.Date(fecha_inicio) + semanas_lag[i] + 8))
plot(estados, add = TRUE)
legend("topright",
c("Número de pruebas por COVID-19 con base en los reportes SSA de México"))
legend("bottomleft",
as.character(pruebas_breaks_labels),
fill = brewer.pal(n = 9, name = 'YlGn')[1:9],
title = "Por 100 mil habitantes"
)
## grafica Pruebas
plot(mis_municipios,
col = mis_municipios$Color_raw,
axes = TRUE,
border = NA,
#col = "transparent",
main = paste0("Pruebas semanales ",
as.Date(fecha_inicio) + semanas_lag[i] + 1," - ",
as.Date(fecha_inicio) + semanas_lag[i] + 8))
plot(estados, add = TRUE)
legend("topright",
c("Número de pruebas por COVID-19 con base en los reportes SSA de México"))
legend("bottomleft",
as.character(pruebas_raw_breaks_labels),
fill = brewer.pal(n = 9, name = 'PuBuGn')[1:9],
title = "Pruebas semanales"
)
# ## grafica defuunciones indice
# plot(mis_municipios,
# col = mis_municipios$Color_indice,
# axes = TRUE,
# border = NA,
# #col = "transparent",
# main = paste0("Índice de desviación del número de pruebas semanales esperado ",
# as.Date(fecha_inicio) + semanas_lag[i] + 1," - ",
# as.Date(fecha_inicio) + semanas_lag[i] + 8))
# plot(estados, add = TRUE)
# legend("topright",
# c("Áreas en blanco reprezentan zonas sin pruebas realizados en el periodo",
# "Número de pruebas COVID-19 con base en los reportes SSA de México"))
# legend("bottomleft",
# as.character(pruebas_indice_breaks_labels),
# fill = brewer.pal(n = 9, name = 'PRGn')[1:9],
# title = "Indice de pruebas"
# )
# ## grafica defuunciones indice
plot(mis_municipios,
col = mis_municipios$Color_positividad,
axes = TRUE,
border = NA,
#col = "transparent",
main = paste0("Positividad en pruebas de COVID-19 ",
as.Date(fecha_inicio) + semanas_lag[i] + 1," - ",
as.Date(fecha_inicio) + semanas_lag[i] + 8))
plot(estados, add = TRUE)
legend("topright",
c("Áreas en blanco reprezentan zonas sin pruebas realizados en el periodo",
"Número de pruebas COVID-19 con base en los reportes SSA de México"))
legend("bottomleft",
as.character(pruebas_positividad_breaks_labels),
fill = brewer.pal(n = 9, name = 'Reds')[1:9],
title = "Positividad"
)
## verification plot wit coordinates only
#plot(x = mis_municipios@data$Longitud,
# y = mis_municipios@data$Latitud,
# col = mis_municipios@data$Color_raw,
# type = "p", pch = 19)
}
knitr::asis_output(htmltools::htmlPreserve(cc_html))
.